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We revisit the issue of the large negative next-to-leading order (NLO) cross section for 
single inclusive hadron production in p A collisions in the saturation formalism. By im¬ 
plementing the exact kinematical constraint in the modified dipole splitting functions, two 
additional positive NLO correction terms are obtained. In the asymptotic large k± limit, we 
analytically show that these two terms become as large as the negative NLO contributions 
found in our previous calculation. Furthermore, the numerical results demonstrate that the 
applicable regime of the saturation formalism can be extended to a larger k± window, where 
the exact matching between the saturation formalism (in the asymptotic k± regime) and the 
collinear factorization calculations will have to be performed separately. In addition, after 
significantly improving the numerical accuracy of the NLO correction, we obtain excellent 
agreement with the LHC and RHIC data for forward hadron productions. 
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I. INTRODUCTION 

A focal point of the frontier of high energy nuclear physics at RHIC and the LHC is the study of 
saturation in hadron collisions. Saturation is an effect that emerges due to bremsstrahlung gluon 
radiation in the hadronic wavefunction. It was prompted by the theoretical prediction that, at 
high energy, the gluon density rises rapidly as x, the longitudinal momentum fraction of the gluons 
with respect to their parent hadron, decreases. This rise is governed by the famous Balitskii, 
Fadin, Kuraev, and Lipatov (BFKL) evolution equation |1], which emerges from resunnnation 
of terms proportional to a s ln^. However, when the gluon density becomes high, it is expected 
that gluons start to recombine and QCD dynamics becomes nonlinear. This eventually leads to 
the onset of gluon saturation M, as a result of the nonlinear QCD dynamics. To quantify the 
recombination effect, a nonlinear term in the gluon evolution equation was proposed in Ref. ms]. 
This nonlinear extension of the BFKL evolution equation was later independently derived by 
Balitsky and Kovchegov; accordingly, the equation is referred to as the BK evolution equation UM- 
Theoretically, it seems that gluon saturation is bound to occur as a result of high energy evolution. 

The task remains to find a “smoking gun” signature of gluon saturation in experimental data at 
e.g. RHIC or the LHC. A wealth of results in p A collisions, ideal for observing saturation [7], are 
becoming available. But it is critical to have precise, quantitative phenomenological calculations 
in the saturation formalism to compare to these experimental results. 

Single inclusive hadron production in p A collisions at high energy reveals the interesting physics 
of gluon saturation particularly well, compared to pp collisions. The effect of dense gluons in the 
target nucleus can be characterised by the introduction of a semi-hard momentum scale, which 

is known as the saturation scale Q s , a function of the momentum fraction x g = k + luon , and the 

"'nucleon 

nuclear mass number A. Roughly speaking, the saturation scale can be used to separate the 
saturated (dense) regime, in which the nonlinear energy evolution applies, from the (dilute) regime 
in which the evoution is linear. When the typical hard scale of the scattering, Q, is less than Q s , 
one expects that the target partons involved in the interaction are saturated. On the other hand, 
when Q Q s , the saturation effect is no longer important, and standard perturbative QCD should 
be sufficient to describe the data. 
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It is generally believed that the transverse momentum of typical gluons inside nuclear targets 
is roughly Q s (x g , A). In high energy pA collisions, before partons from the proton projectile 
fragment into hadrons with transverse momentum p ±, they undergo multiple interactions with the 
dense gluonic fields in the highly boosted target nucleus, picking up a transverse momentum of 
roughly Q s in the process. The squared saturation momentum Q ^ in nuclear targets is A 1 / 3 times 
of that in protons, due to random multiple scatterings. Therefore, the transverse momentum (p±) 
spectrum of the produced hadrons exhibits different behaviour in p A collisions, especially in the 
relatively low p± regime, than in pp collisions. 

Measurements of single inclusive hadron production in pA collisions at RHIC j8hT3] and the 
LHC |14H21j have provided plenty of data. There have been many theoretical and phenomenological 
efforts [22fl39j on this subject, which imply that gluon saturation (or shadowing) plays an important 
role in the production of forward rapidity hadrons. 

The first complete next-to-leading order (NLO) calculation for inclusive hadron productions ffO] 
SB in p A collisions was achieved a few years ago, using the Mueller’s dipole formalism (see Ref [22j 
for discussion at NLO). This computation is the first complete NLO calculation which uses di¬ 
mensional regularization and with the MS regularization scheme, which is necessary to correctly 
incorporate the available NLO collinear parton distributions (PDFs) and fragmentation functions 
(FFs) without introducing additional scheme dependence. 

Schematically, the full NLO cross section for hadron production at forward rapidity y, with the 
hadron having transverse momentum p± = zk± , can be written as follows: 

dy<¥p_L = J x pf i ( x p}® Dh /^ z }®Fi 9 {k 1 _)®U {0) + ^ J x p fi(x p )®D h/j (z)< ^ , (1) 

All the hard factors T~L are given in Refs. m 23 ! • ^(k ± ) and which are functions of the 

transverse momentum k± of the produced parton, represent the Fourier transforms of dipole scat¬ 
tering amplitudes. x p fi(x p ) and D h /j(z) are the parton distribution and fragmentation function, 
respectively. 

The first numerical analysis of forward hadron production in p A and dA collisions in the small-x 
saturation formalism, based on the NLO results in Refs. mm, was presented in Refs. [HSjEj. 
It was found that the theoretical uncertainty is significantly reduced compared to leading order 
(LO) results, and the calculated NLO cross section agrees well with forward-rapidity RHIC data 
for p_l < Q s - Recall that Q s is the characteristic scale for the gluon density in a heavy nucleus. In 
general, the p± region in which the calculation and results agree increases with the center-of-mass 
scattering energy y / s/v)v, since the typical gluon density probed is larger at high energy. However, 
the numerical results of the NLO cross section abruptly drop to negative values above some cutoff 
momentum which is generally slightly greater than QsQ 

Strictly speaking, the saturation formalism always takes the high energy limit s — > oo, which 
yields large saturation momentum Q s . In this limit, the NLO results in Refs. mm are obtained 
after the subtraction of the rapidity divergence (associated with the small-x evolution in the s —> oo 
limit) as well as the collinear divergences (associated with the DGLAP evolution of PDFs and 
FFs). However, in phenomenological studies of saturation physics, the center-of-mass energy of 
scatterings is finite and the saturation momentum is not very large. It is natural to expect that 
the saturation formalism works when k± < Q s . On the other hand, when k± ;§> Q s , the saturation 
formalism is believed to be no longer applicable since x g is no longer small, and the collinear 
factorization approach should be the relevant formalism to describe the large k± part of the cross 

1 The relationship between the saturation scale Q s , the cutoff momentum at which the results become negative, 
and the boundary of the region in which the calculation accurately describes the data appears to be some sort 
of rapidity-dependent proportionality, but the details are not clear. Ref. [44] includes some discussion of the 
relationship among these momenta. 
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section. Studies [45] have shown that the transverse momentum dependent (TMD) factorization 
is closely related or equivalent to the small- 2 ; factorization in terms of gauge links [461448] . if one 
puts them in the same kinematical region. A similar kinematic restriction is also required for the 
TMD factorization in the hard scattering processes, where the transverse momentum k± is much 
smaller than the hard momentum scale Q (like the invariant mass of lepton pair in the Drell-Yan 
process) [59]. In other words, the TMD factorization would be invalid in the large k _~ Q region. 
In this region, a matching to the collinear factorization calculation is usually performed [49]. This 
argument applies to the case studied in this paper as well. 

In Ref. [ 50] . the matching between the perturbative results is reached by taking the large k± limit 
of the NLO cross section in the small-x formalism and enforcing the exact kinematic constraint. 
One might wonder if there is a way to naturally implement the kinematic constraint in the small-x 
formalism. Doing so could help extend the applicability of the formalism to the large k± regime. 
Recently, the authors of Refs. [5DI5U performed another independent NLO calculation of single 
inclusive hadron production in pA collisions. Their discussion of the Ioffe time dependence, which 
is equivalent to the kinematical constraint, motivated us to investigate the details of the dipole 
formalism application in this process and compute the effect of the kinematical constraints. We 
find that we obtain two additional NLO corrections from incorporating these constraints. These 
two terms were conjectured to be small in the s —> 00 high energy limit, and therefore implicitly 
neglected in the original derivation of the NLO corrections mi im. It is interesting to note that, 
in Ref. [53], a similar logarithmic term played an important role in deriving the Sudakov factor in 
other hard processes. We need to emphasize that there is no Sudakov factor in the process of single 
hadron productions. As we will show later, these additional NLO corrections are indeed small as 
long as k_ j_ < Q s . However, they become large when k± > Q s ■ By including these two terms, we 
can offset the negativity of the NLO terms at larger p± and extend the applicability window of the 
saturation formalism towards larger p± for single hadron productions. 

The goal of this paper is to revisit the issue of the negative NLO cross section found in the large 
k _l regime of forward-rapidity single inclusive hadron production in p A collisions, and the numerical 
implementation of the exact kinematics in the small-x formalism at one-loop order. We find that, 
with the exact kinematical constraint imposed, we can obtain two additional NLO hard factors, 
with one from the q —>• q channel and the other from the g —>• g channel. We first analytically show 
that these two terms are large enough to partially overcome the negative NLO terms found earlier 
in the simple Golec-Biernat and Wusthoff (GBW) model [ 54] , Numerically, these two terms are 
found to be negligible when k± < Q s , but they become important when k± rises to ~ 2 Q s and 
higher. 

More importantly, we have significantly improved our numerical implementation of all NLO 
corrections, which allows us to do phenomenological studies at the LHC energy. We find excellent 
agreement between the full NLO cross section and the forward hadron production data at the LHC. 
This paves the road for a quantitative and precise phenomenological test of saturation physics at 
the LHC. 

The rest of this paper is organized as follows. In Sec. II, we present a detailed derivation of 
the implementation of the kinematical constraint into the dipole model, and obtain two additional 
NLO corrections (one for the quark channel and one for the gluon channel) after subtracting the 
corresponding small-x large logarithms. We further evaluate these two terms by Fourier transform, 
and demonstrate that they have the same behavior exhibited by perturbative QCD in the high 

k ± 

k± limit. In particular, we analytically show that the additional terms are large compared to the 
negative cross section at NLO. In Sec. Ill, we present the numerical results and compare them to 
RHIC and LHC results. We summarize our paper in Sec. IV. 
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II. DIPOLE MODEL AND KINEMATICAL CONSTRAINTS AT NLO 


Let us first review the exact kinematical constraint discussed in previous publications 
since it plays an important role in this paper. We will continue to use light-cone perturbation 
theory as in Refs. B00H, and define p + = - JT and p = - ^ . 

The kinematical constraint is derived from the conservation of the — component of the four 
momentum before and after interactions for 2 —> 2 processes. For quark production with transverse 
momentum k±, as illustrated in Fig. [lj we obtain 

l 2 k 2 

+ ( 2 ) 


2(1 -£)x p P+ 


2 £x p P+ 


As £ approaches 1, the above kinematical constraint indicates 


l± < (1 ~OxpS, 


or 


£ < 1 -— ■ 


x n s 


( 3 ) 


For a rapidity divergent term, we find that the above constraint modifies the upper limit of the 
divergent integral as follows 


f 

■Jo 


1 x p s 1 k ± 

- -- = In ~2~ = In — + In - jT , 

1 - (, l \ Xg Z | 


( 4 ) 


where x g = according to the leading order kinematics. In the high energy limit, we assume 

k _l ~ Zj_ which makes the second term very small. However, as we shall show in the following 

k 2 . 

discussion, the In term becomes important when k± gets larger than the typical saturation 

( _L 

k 2 

momentum. The direct evaluation of this In term is not easy in the momentum space, but 

l ± 

indirectly evaluating it in coordinate space is quite straightforward, as we shall demonstrate in the 
following discussion, after encoding it into the modified dipole splitting functions^ 

Inspired by Refs. EHEZI, we find that it is convenient to encode the kinematical constraint by 
modifying the dipole splitting function for the q —>• qg splitting, shown in Fig. [lj as follows: 


^(p + ,k + , U± ) = 2ttL 


(i - 0 p + 


[1 - J 0 (u±A)] < 


u±-e 


( 1 ) 


u±-e 


'( 2 ) 


-(5 a -8/3- + £5 a _|_5g-|_), A — 1, 
■(^a+^/3+ + £8a-8p-)i A = 2. 


( 5 ) 


2 A similar term, In , where M is another hard scale (e.g. the Higgs mass or dijet invariant mass), is used to 
L ± _ 

derive the Sudakov factor when M 2 k\ [S3]. The Sudakov physics is different from what we are discussing here, 
since basically the single hadron production process is a single scale problem. First of all, the new terms that 
we obtain can be never interpreted as a Sudakov factor, instead it should be viewed as part of the NLO power 
correction. In addition, the color flow is completely different. Note that the color factor for the quark Sudakov 
factor is CV while the one for BFKL physics is always N c / 2. 
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with A 2 = £(1 — £)x p s. The original dipole splitting function, which is proportional to 


( 1 , 2 ) 


arises from the Fourier transform of 


h 


A 1,2) 


The additional term of — Jo(u±.A) arises from the 


kinematical constraint, Eq. ([3]), imposed during the Fourier transform. 


In general, the correction to the splitting function Jo (u±A) does not play any important role 
since it vanishes when we take the high energy limit s —> oo. It only becomes important when the 
gluon longitudinal momentum fraction 1 — £ approaches zero. Specifically, for quark production at 
one-loop order, we always get the following DGLAP-type contribution from real diagrams: 




1 + £ 2 

( 1 - 0 + 



( 6 ) 


For the first term on the right hand side of the above equation, we can safely take the s —)> oo limit, 
since this term is regular when £ -+ 1. However, one can not neglect the correction Jo ( u±A ) for 
the second term, since it is singular when £ —> 1. Clearly, in this NLO calculation for single hadron 
production in pH collisions, the kinematical constraint only affects the rapidity subtraction term. 

The relevant contribution to the cross section, with the modified splitting function, can be 
written as 


a 


f 1 d £ f d 2 x±d 2 y ± d 2 b 1 _ _ ikl . (xi - yi ) 


2vr 2 


/o i-e 

X 


(2vr) 2 

l2 


—S(x±,y±) + S(x±,bj_)S(b±,y±) 


f [!-»+)] + [1-W] _ ^ (1 _ Jo( „ ±A)] (1 _ JoK A,]}, 

( u ± ur u.uf J 


(7) 


where u± = x± — b± and u' ± = y± — b± J^] One can actually approximately integrate over £ and find 
that 


d£ 

i-£ 


1 - J 0 (u_L \Jx p s{ 1 - £)) 


In 


XpSUj 


= In-b In 


k 2 , u 2 , 


x n 


f 1 d£ 

/o 1-^ 


1 - J 0 (u±^jx p s( 1 - £)) 1 - Jo^u^XpS^ - £)) 

, x„su±u\ ,1 k 2 , u±u 'I 

~ m ——^—— = In-b In 9 , 


( 8 ) 

(9) 


with co = 2e~ lE . Here we have used x p x g s = k\. It is then clear that the first term In y- can be 
subtracted from the NLO cross section and interpreted as the BK evolution of the dipole amplitude 
S up to rapidity Y g = In d-. The second term in the above equations, which is conjugate to the term 

In jA as in Eq. Q (see also Eq. (3.12) of Ref. [5.2]) with l± being the gluon transverse momentum, 
arises due to the exact kinematical constraint. More precisely, the second integral should give 
In — + In instead of In — + In k l u ^ a ±- w ith u < = min 

Xg Cq Xg Cq 

for L q (k±) non-analytical and the precise numerical evaluation more challenging. Fortunately, one 
can numerically check that the resulting L q (k±) has similar large-/c_|_ behaviour, and it gives the 
same high k± tail, since u± ~ u' ± when k± —> oo. Besides, as we will show later, in the low-/cj_ 
region, L q {k_ |_) is negligible in the total cross section. Our goal here is to extract the correct large 
k± tail of the additional hard factor, which eventually helps to extend the applicability of the 


{ri_L,which makes the calculation 


3 This term looks similar to the last term of Eq. (4.21) in Ref. [52]. Here we are implementing the kinematical 
constraint, while they are discussing the Ioffe time cutoff. These two become equivalent if one identifies their 
as the center-of-mass energy s in our paper. 
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small-x calculation. In this sense, we can approximate In 


± u '< 

s 

-o 


as In 


k 2 ?/. | u f 
c 0 


Also, because the 


rest of the expression is symmetric under the exchange u± <->• u\ , it is useful to note that this is 
equivalent to using In 


~P r 

c o 


or In 


The leftover terms can be cast into an additional hard factor which reads 


L q (k ± ) = 


a s N c 

2i r 2 


d 2 x_i_d 2 y_i_d 2 6_i_ 
—e 


ik 


(27T) 


[ 5 ( x± - b ± )S(y± - b ± ) - S(x± - y ± )] 


1 , kf uf 1 , k 2 u'f 2u± - 7,2 


—o- In 
u\ 


+ ~i2 In' 

uf 


uf u'f 


In 


kf\u±\ 


u , 


( 10 ) 


The corresponding contribution to the single inclusive cross section in this channel can be written 
as 


d 3 a Lq 

dyd 2 p± 



( 11 ) 


It is not hard to show that the above contribution from L(k±) is free of both UV and IR divergences. 
When 6 _l —> x±, the first bracket vanishes. When b± —> 00 , the second bracket vanishes. Due to 
these strong cancellations, it was believed that this contribution should be small. In fact, Ewerz 
et al [55] studied the Ioffe time effect of the dipole model in deep inelastic scattering for inclusive 
total cross sections, and they found that this effect is small. For single inclusive hadron production 
in p A collisions, as we demonstrate below, the effect is small when p± is small, but it becomes as 
large as other NLO corrections when p± ~ Q s . 

Note that this term is physically and fundamentally different from the so-called AH correction 

f,2 

from Kang et al [56] , which is proportional to the rapidity interval Y — Y g = In — + In . As 
commented in Ref. the choice of the rapidity interval leads to an unphysical conclusion and 
violates the small-x factorization. The new additional term L q (k±) does not depend on either the 
projectile longitudinal momentum fraction x p , or the hadronic mass m p . It is important to notice 
that QCD factorization does not allow us to have hadronic mass rn p in any hard factors. Otherwise, 
this implies that we can not separate the non-perturbative physics from the perturbative calculable 
hard factors. It is also clear from our above derivation that x p naturally cancels out and thus does 
not appear in L q . We would like to emphasize that the so-called AH correction discussed in 
Ref. [56] is unjustified and should be absent in view of the small-x factorization. 

Let us derive the following simplified expression for L q (k±) which is much easier to evaluate 
numerically. It is straightforward to use the following Fourier transform identified 




ln k ± u ± 
2 111 2 

U 1 c 0 


1 

87T 


d 2 l±e U ± - U± 



1 

2tt 


d 2 l ± e il± ' u± 


il _l 


ln 


k± 

ll 


( 12 ) 

(13) 


4 Note that this Fourier transform may be problematic for u± = 0, therefore we should exclude the point where 
u± = 0, in principle. However, since the first bracket in Eq. ( |10[ ) vanishes when x± —¥ b± (or equivalently 
ui —> 0), which suggests that we are justified in ignoring the fact that L q (k±) is undefined at that point. We have 
also numerically tested that the two expressions of L(k±), before and after the Fourier transform, give the same 
numerical results. 
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to find that 


a, N r 


Lq(k±) — ^_2 '-’-L [Lql(k±) + L q2 (k±) + L q3 (k±)\ 

- ik ^S(r ± )( In^p) 2 , 


with 


Lql(k±) = - 


d J± e -i 

2n 


Lq2 (k ± ) = (2vr )F(k ± ) j d 2 l ± F(k ± - l ± ) (in 


L q3 (k ± ) = -4 / d 2 /_i_d 2 / / _ L F(/cj_ - 4)F(fe ± - Z ± ) 




l± In k * 


m 2 l 2 


(14) 

(15) 

(16) 
(17) 


In deriving the above expression, we have used the fact that the impact parameter integration 
simply gives S±, which is the area of the target nucleus, (see Appendix for detailed derivations.) 
In fact, one can further evaluate L q (k±) analytically in the GBW model by assuming 


S(R±) = exp ( — 


QlR 2 


F(k ± )= 


l-R 


ttQ 


2 exp 




Qs 


(18) 


and hnd 
L 


,i(U) = ^ 1 


— 1 — 

’ Qi 


-2 In 


k‘j_e YE 

~qT 


L (b°) 


k 2 1 

-i 

’ Qi 


+ 


In 


k 2 , e lE 


L q2 (k±) = dn 2 F(k±)F(k±) / d/j_/j_exp 


L q s(k x ) = ~F(k ± ) 


l - eXp[ -Ql 


*OC 

dl±h 


Qi 

2/j_/c_i_ 

~qT 

21 ± k 

~Ql 


2 7T 2 ' 

+ Y 


exp 



, k i 

If 


exp 


/ 2 

Qi 


, k l 

ln IT 


(19) 


For L q2 (k± ) and L q3 (k±), in principle, one can also perform the dZj_ integration and obtain ana¬ 
lytical final results in terms of derivatives of hyper-geonretrical functions. Asymptotically, L q \{k±) 

and L q3 (k±) give -pr and — ~pr in the large k± limit, respectively, while L q2 (k± ) is exponentially 


suppressed. The most interesting observation is that L q (k_ j_)| 
the NLO quark channel hard factor, which involves 

asNc f d ,^ (1 + ^ asNc 

't/z 


^ a s N c S± 4 Q 2 S 
k±—>oo 47r 2 A:'j 


47T 2 


(1 - 0+ ki 


X q( x )^1 9l 

4vr 2 XpQ[Xp) 12 k 4 ! 


Comparing to 


( 20 ) 


in the large k± limit, it is conceivable that this term is sufficient to largely cancel the large and 
negative NLO cross section found earlier. In order to understand the importance of the additional 

contribution as in e T @> it is illuminating to compare it with the leading order 

section in the quark channel which can be written as 


cross 


d 3 


a 


LO 


dyd 2 p ± 


= J z 5 ^ 2 xpqf(xp)D h/q (z)S ± F(k ± ), with F(k±) = j lk± ' r± S(r ± ), (21) 


since the LO cross section can provide an order of magnitude estimate of the total cross section. 
It is then straightforward to see that one just needs to compare L q (k±) (after factorizing out S ±) 
with F{k |_). As shown in Fig [2j L(kj_) is small compared to F(kj_) in the low kj_ region, therefore 
the additional contribution is negligible when k± < Q s - One the other hand, L q (k± ) falls slowly 
with k± and becomes important when k± > 2 Q s . 
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• 10" 4 




k_L /' Q s 


kj_ IQs 


O 2 

FIG. 2. The comparison between ^ L q (k ±) and QgF(k±) with Sj_ factored out. (One can also simply set 
Q s = 1 GeV.) Here we have employed three different numerical methods to evaluate L q (k±). The blue dots 
indicate the direct numerical evaluation of L q (k±) as in Eq. (101 with Mathematica, while the red circles 
represent the evaluation of Eq. (14). The golden diamonds correspond to the numerical results obtain from 
Eq. (141 by using our SOLO code programmed with C + +. The asymptotic k± behaviour of L q (k±) is 
indicated by the green dashed line. The numerical uncertainties are very small as shown in the right plot. 


III. THE GLUON CHANNEL CONTRIBUTION 


For the gluon channel, we can use the same procedure to take into account the kinematical 
constraint. As discussed above, all the hard factors were computed in Refs. mum except for the 
rapidity divergent piece. For a large nucleus target and in the large N c limit, we can write the 
rapidity divergent term as 


a,N, 


7r 


Sl,c d2X tf b± e ~ ik±<X± ~ V±) [~S(x±, y+) + S(x ± , b ± )S(b x , y ± j] S{x ± , y ± ) 


x 


[1 - Jo (uj_A)] 2 [1-J 0 «A)] 2 2u ± -u' ± 


+ 


w 


u 


1 2 


u 2 , u'? 


[1 - Jo («LA)] [1 - Jo KiA)] \ , (22) 


where u± = x± — b± and v! j_ = y± — b±. We have put in the kinematical constraint in the 
same fashion as we did for the quark channel. Again, after subtracting the small-a: logarithm In A- 
through the BK evolution equation for the dipole scattering amplitude in the adjoint representation 
appearing in the gluon channel, we can obtain the remaining additional hard correction due to the 
kinematical constraint 


L g(k+) = °^ [ [-S(x ± ,y x ) + S(x ± ,b ± )S(b ± ,y ± )]S(x_ L ,y_ L ) 


7T 


—q- In 
IF, 


k 2 ,u 2 , 1 , kill'? 2u±-u\ kl |it / i I 

m - 


+ ~12 ln 


u 


w 2 u’ 2 


(23) 


At the end of the day, we can find the following additional contributions from the gluon channel 


dyd 2 p± 


d z 

^ x p9(xp)D h /g{z)L g (k_ l), 

r % 


(24) 
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with L g (k ±) = g ^ L S _l [L 5 i(A:_l) + L g2 {k±) + L g3 (k ±)] and 

^gi(fc±) = -2 J < ^^e~ lk± ' r± S(r ± )S{r±)(^n l ^p^ S j , (25) 

i ff2 (fcj_) = (4 tt) J d 2 l ± d 2 l' ± F(k ± -l±- l'±)F(k± - 4)F(4) (in ^ , (26) 

£ s3 (fc±) = -8 J d 2 l L d 2 l’ L d 2 l" L F{l'[)F{k L -l' ± - l'i)F(k ± -l ± - In ^ (27) 


In the GBW model, using Eq. (18), these expressions can be simplified to 


L„1 (*±)= -Qj( i(2,0) 


- 1 ,-: 


2Q2 




205 

+ 


- 1 ,-: 


K 

2QI 


In 


k 2 t e lE 


2 Q 


2 vr 2 ' 

+ y 


exp 


k 2 


2 Q 


167T f°° 

Lg2{k±) = - ^F(k±) / dl ± l ± / dq±q ± I 0 

Vs J o Jo 


2k ± q ± \ T f2l ± q ± \ -li+^1 ^ ^.2 

Q2 JH 02 J e * Pi 


(28) 


(29) 


L g3 (k ± ) =-^F(k ± ) dl ± J^ dq ± (l-e-^ Q ^I 0 (^^yi(^^y In 


k i 

l 2 , 


(30) 


It is straightforward to find that the leading power behaviour of L g (k±) is a °F, c S± with 
defined as the quark saturation momentum. 


IV. NUMERICAL RESULTS AND PHYSICAL DISCUSSIONS 

A. Numerical setup 

The calculation of single inclusive hadron production in pA collisions up to NLO level with the 
full running coupling has been recently implemented as the computer program Saturation physics 
at One Loop Order, or SOLO [43]. Initial results from the program showed that the full NLO 
pA —> hX cross section agrees fairly well with the forward RHIC data. The same paper also 
provided numerical results at the energy scale of the LHC at extreme forward rapidities y > 5, 
since at the time no LHC experiment had published any results. 

Since then, ALICE, CMS, and ATLAS have released their experimental data of the pA —> hX 
differential cross section, all of which are at roughly central rapidity, —2 < y < 2. Unfortunately, 
the initial version of SOLO gives results with very large uncertainties for LHC central rapidity 
collisions, for basically the following reasons: first, several of the terms involve oscillatory factors 
of the form Jo( Px J x ) ° r Jo( Px g X ); which are integrated over r±. At the LHC, the physical scenario 
where a high x p parton from the proton projectile radiates a soft (small £) gluon (or quark), which 
then fragments into the produced hadron, becomes a much more significant contribution than at 
RHIC. Specifically, £ or z can be as small as r = , which decreases from 0.04 at BRAHMS 

with y = 2.2, to 0.0002 at ALICE or ATLAS with y = 0 for p± = 1 GeV. Such small values of 
z and £ produce rapid oscillations in the integrand, which most numerical integration algorithms 
are notoriously bad at handling. Although there are specialized algorithms available, they are 
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prohibitively difficult when the integrals have 4, 6, or even 8 dimensions, as is the case with most 
terms in the cross section. 

Other terms involve factors like F(k± — 1±) which are integrated over l±, or similar factors with 
other forms of the argument. These functions have their peaks where the argument is zero: for 
example, l± = k± in the first case. The numerical integration algorithms used by SOLO are most 
effective when the function being integrated has its peak where the integration variable l± = 0, 
which is approximately satisfied at RHIC (the largest contributions come from k± ~ lOGeV), but 
not at the LHC. Under LHC kinematical conditions, we have to accommodate much larger values 
of k _l while keeping the numerical uncertainties under control. 

Through a series of transformations and various tricks, together with a tremendous amount of 
effort, we have converted the formulas originally used by SOLO |44| to improved versions with vastly 
smaller numerical uncertainties. At this moment, we believe that all the numerical uncertainties 
are well under control even at the middle rapidity region of the LHC. Appendix B explains some 
of the key points of these transformations, but we will reserve the full details of the numerical 
implementation for a separate document to be released later. 

We have computed results for two different parametrizations of the dipole scattering amplitude 
S (or its momentum space expression, F). First, the GBW model, defined by 

Sgbw(cl) = exp^-^j^j (31) 

with 

Qs(x g ) = cA l/z Ql (^j (32) 


where A represents the number of nucleons in the target nucleus and c = \J\ — b 2 /R 2 , where b is 
the impact parameter and R is the nuclear radius, accounts for the centrality of the collision. We 
use the same parameters as in previous SOLO results [4311441150| : c = 0.56 to represent minimum 
bias collisions, and the original GBW fit parameters from Ref. [54], xq = 0.000304, A = 0.288 and 
Qq = 1 GeV 2 , which are based on HERA data. The GBW model is often used in phenomenological 
calculations due to its simple analytical form. 

In addition, we also show results computed using the numerical solution of the BK equation 
with the running coupling correction j58l - lfil| . setting the QCD running coupling scale in the rcBK 
equation to A = 0.1 GeV. This solution has been shown }62H65| to be very useful in phenomenol¬ 
ogy, especially at high transverse momenta. Previous studies have also considered the McLerran- 
Venugopalan (MV) model |66j . an analytical expression which gives a high-/cj_ power tail similar 
to the BK solution, but we have omitted those results from the present analysis because the GBW 
and rcBK models are sufficient to show the interesting features of the results. 

Results computed from SOLO are for the quantity ULaFTV”’ can trivially converted 

into the differential yield for production of a single pion species in the center-of-mass frame, 

d 3 N pA^x _ s± 1 d 3 a pA^x 

dyd 2 p± cr ine i S_ L dyd 2 p± 

where n is a single species of pion: tt + , 7r°, or tt~ . However, the experiments measure 


1 d 2 N pA ^ hX 1 d 2 a pA ~^ hX 

2np_L dr/dp± 27rp_Lcr ine i d?/dp_L 


(34) 


where h may include several different hadron species, depending on the detector, and cxi ne i is the 
total inelastic cross section. We neglect the difference between rapidity in the lab frame and 









11 


y = o y = 0 

(RHIC) (LHC) 


y < o 


y> o 


p,d 


A 



forward 

hadron 

production 


FIG. 3. The orientation of rapidities used by SOLO and throughout this paper. Positive (or forward) 
rapidity is always in the direction of the proton (or deuteron, for RHIC) beam. Some results published by 
ALICE [T5] and ATLAS 21] use the opposite orientation. In this paper we always use y to represent the 
rapidity in the center-of-mass frame. 


pseudorapidity rj. Accordingly, we multiply the output from SOLO by a factor ^S±/a me \ to make 
it compatible with the experimental measurements. 

• BRAHMS measures negatively charged hadrons. We use ^ = 1.3 to account for the yields 
from kaons and other hadrons. We also set S± ~ 7r(7.5fm) 2 = 1770 mb, and use cTi ne i = 
2400mb [67j. 

• STAR measures only neutral pions, so ^ = 1, and <7i ne i = 2210 mb mi. with the same S± 
as BRAHMS. 

• ALICE and ATLAS measure all charged pions, kaons, and protons. We use a result from 
CMS [68] that the kaon and proton yields are 13% and 6%, respectively, of the pion yield, 
giving ^ « 2.4, and cxi ne i = 2100 mb from LHCb [69]. For lead nuclei, S± « S^ u x 
(208/197) 2 / 3 = 1830 mb. 

As far as the definition of rapidities is concerned, deuteron beams have positive rapidity and gold 
nuclei beams have negative rapidity at RHIC. The energy of both beams per nucleon is y/~SNN = 
200 GeV. Therefore, the center-of-mass frame is the same as the lab frame. The specification 
of rapidities in the SOLO package follows the above setup: positive rapidity always refers to 
the deuteron-going direction (or proton-going, at the LHC), as shown in Fig. [3j On the other 
hand, some of the ATLAS and ALICE pPb data measured at y/sjyN = 5.02 TeV are presented 
in the opposite rapidity configuration, with the proton beams having negative rapidity. In order 
to compare our results with those data without confusion, we have flipped the sign of rapidity in 
the experimental results and labeled our plots with the rapidity y in the center-of-mass frame in 
the SOLO convention. Throughout this paper, forward rapidity, y > 0, always means the rapidity 
region along the proton (or deuteron) beam direction in the center-of-mass frame. 


B. Discussion of numerical results 

Figures [4] and [5] show the differential dA —> hX yields at forward rapidity for BRAHMS [9j 
and STAR [10], respectively, along with the corresponding results from SOLO using the new 
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P±[ GeV] 


P±[ GeV] 


FIG. 4. Comparisons of BRAHMS data [9] with the center-of-mass energy of v / s]\nv = 200 GeV per nucleon 
at rapidity y = 2.2,3.2 with our results. As illustrated above, the crosshatch fill shows LO results, the 
grid fill indicates LO+NLO results, and the solid fill corresponds to our new results which include the NLO 
corrections from L q and L g due to the kinematical constraint. The error band is obtained by changing /i 2 
from 10GeV 2 to 50GeV 2 . 


(transformed) formulas. The LO and LO+NLO curves are very similar to earlier results published 
in Ref. [l3|; some slight differences are due to the increased precision of the new formulas. In the 
meantime, the L q and L g corrections are completely negligible in the region where p± < Q s . On 
the other hand, where p± > Q s , L q and L g start to become important and alleviate the negativity 
problem in the GBW model, and help us to better describe the data in the high p± region. In the 
rcBK case, we find that the full NLO cross section now becomes completely positive and provides 
us excellent agreement with all the RHIC data. 

In Figure [6j we show the comparison between the forward ATLAS data at y = 1.75 and the 
numerical results from SOLO. We observe remarkable agreement between the full NLO calculation 
from the saturation formalism and experimental data up to 6 GeV. Again, as we have seen earlier, 
the newly added L q and L g corrections help to increase the applicable p± window of the saturation 
formalism from roughly 2.5-3 GeV to 6 GeV. From 6 GeV and up, the full NLO cross section 
still becomes negative, which implies that the saturation formalism does not apply anymore and 
the collinear factorization should be used. Admittedly, what we have seen is only one piece of 
a promising clue for the gluon saturation phenomenon. More data in different forward rapidity 
windows at the LHC would allow us to conduct precise tests of the theoretical calculation, and 
may eventually provide us the smoking gun proof. 
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FIG. 5. Comparison of STAR data JO] with snn = 200 GeV at y = 4 with results from SOLO for the 
GBW and rcBK models. The color scheme is the same as in figure [4j and again, the error band comes from 
fj, 2 = 10 GeV 2 and 50 GeV 2 . We do not see the negative total cross section because the cutoff momentum 
above which the cross section becomes negative is larger than the p± of the available data, and in fact larger 
than the kinematic limit ^/s^nS-~ v ■ 



fv[GeV] P_L[GeV] 

FIG. 6. Comparison of ATLAS forward-rapidity data m with the center-of-mass energy of y'sjvjv = 
5.02 TeV at y = 1.75 with SOLO results for the GBW and rcBK models. Again, the color scheme is the 
same as in figure [4] Here the error band shows plots for p 2 = 10 GeV 2 and p 2 = 100 GeV 2 . Since the 
numerical data for these measurements are not published, we have extracted the ATLAS points from Fig. 6 
of Ref. [21!. The extraction procedure introduces uncertainties comparable to the size of the points. 


In Figure [7J we show the comparison between the ALICE and ATLAS data at y = 0 and the 
numerical results from SOLO. We find that the full NLO results, especially the one with the rcBK 
solution, miss the data. (It seems that the GBW model roughly agrees with the data, but we believe 
that it is probably just a coincidence.) This indicates that the dilute-dense factorization breaks 
down at y = 0. This is completely expected for the following reason. First, the collinear parton 
distributions of the proton projectile do not resum small-x logarthms and may have considerable 
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rcBK Aq CD = 0.01 



P±[ GeV] P±[GeV] 

FIG. 7. Comparison of mid-rapidity data at y/spfN = 5.02 TeV at y = 0 from ALICE [15j and ATLAS m 
with SOLO results for the GBW and rcBK models. Both data sets correspond to the same kinematic 
parameters. Again, the color scheme is the same as in figure [4] These results display the breakdown of the 
dilute-dense factorization approach when the separation between x p and x g is not sufficiently large. — is 
roughly on the order of 1. 


uncertainties in the very low-x region. Most importantly, the dilute-dense factorization derived in 
Refs. TOl HI] assumes that the proton projectile is dilute while the nuclear target is dense. In the 
forward rapidity region, for example y = 1.75, one can estimate that roughly x g /x p ~ 1CU 2 , which 
indicate that the small-x evolution is more important in the nuclear target than in the proton 
projectile. Therefore, we can use the integrated parton distributions in the proton projectile and 
use the small-x evolution for the unintegrated gluon distribution in the nuclear target. On the other 
hand, in the middle rapidity region, x g /x p ~ 1, we should use the small-x evolution to resum both 
a s In A- and a s In A- simultaneously, since they are of the same order. This means that we need to do 
a complete NLO calculation in the framework of the so-called k± factorization with unintegrated 
parton distributions for both the proton projectile and the nuclear target. Unfortunately, this 
calculation is very challenging, therefore we shall leave it to future studies. 

Some more discussion and comments are in order, as follows. First of all, as we have shown 
analytically and numerically in the GBW model, the two additional NLO corrections derived in 
this work extend the applicability of the saturation physics calculation further into the large p± 
region, without significantly modifying the low p± results. The numerical results using the rcBK 
solution also support the same conclusion. In the very large p± region where the saturation effect is 
extremely small, the full NLO cross section may still become negative, but this is already beyond 
the applicability of the saturation formalism. In this region, it is well-known that the collinear 
factorization is the relevant framework and provides the best description of the QCD dynamics. 

Second, the comparison between the data and our calculation suggests that the implementation 
of the kinematical constraint works slightly better for the rcBK approach. As compared to the 
GBW model, the rcBK solution of the dipole scattering amplitude has the correct perturbative 
tail. 

Last but not least, as shown in the numerical results, the negative cross section at NLO may 
still persist when p± of the measured hadron is much larger than the saturation momentum Q s . 
We recall the results of Ref. |50j: first, that the perturbative QCD calculation from collinear 
factorization can describe data in the large p± region; and also, with exact kinematics (which simply 
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removes the plus function), the perturbative QCD calculation analytically matches the large-pj_ 
expansion of the results from small -x factorization at next-to-leading order in a s . Therefore, the 
nuclear modification factor R p a measured in p A collisions at large p± ( p± Q S (A) > Q s (p )) 
becomes 


= daP A /dyd 2 p ± S a Q 2 (A) = xG A (x) 

N co \\dcrPP/dyd 2 p± AS±Q 2 (p) AxG p (x ) 


Here the number of binary collisions N co n is A in pH collisions. (iV co ii is H 1 / -3 if R p a is computed 
from the measured yields.) This has been repeatedly seen in a lot of experimental data, for 
example PHENIX 0, BRAHMS |9], STAR HD], ALICE pH], CMS HU and ATLAS |3T] . As 
indicated by both the experimental data and our NLO analysis, it seems more and more clear that, 
at sufficiently high energy and forward rapidity, the saturation effect is dominant in the low p± 
region where R p a < 1- On the other hand, the moderate and large p± region, where R p a > 1, 
is described by simple perturbative QCD, with the saturation effects encoded in the subleading 

q2 _ 

power corrections as shown in e.g. Eqs. (28) and (29) of Ref. [70]. To search for clear and 

p± 

compelling evidence of gluon saturation in single inclusive hadron production, one should focus 
on the low-pj_ part of spectra in the forward rapidity region of the proton beam in pA collisions, 
which is dominated by semi-hard scattering in the vicinity of the saturation scale Q s . Having said 
that, one should also be aware that p± of the measured hadron should be kept sufficiently large 
(at least 0.5 GeV) to avoid nonperturbative QCD effects. 


V. CONCLUSION 

In this paper, we have investigated the details of applying the dipole formalism to inclusive 
hadron production in forward pH collisions at the energy ranges of RHIC and the LHC. In partic¬ 
ular, we derived two additional terms by considering the kinematical constraint ([2]) in the dipole 
formalism at next-to-leading order. These two terms were assumed to be negligible in the high 
energy limit s —> oo in previous studies. In order to do more precise and reliable numerical cal¬ 
culations for phenomenological studies of saturation physics, we have to include these additional 
terms at finite center-of-mass energy y/sNjsr- From an extensive phenomenological study, we found 
that these additional terms extend the applicability of the NLO cross section in the small-x satu¬ 
ration physics formalism to higher values of the produced hadron momentum p± . Matching to the 
perturbative collinear factorization will further extend the kinematic coverage of the theoretical 
predictions for this process [50]. 

From the explicit NLO analysis of the single inclusive hadron spectrum in pH collisions, we 
argue that the nuclear modification factor R p a shall approach 1 at sufficient large transverse 
momentum, where the collinear factorization calculations apply to both pp and pH collisions. In 
the low transverse momentum region p± < Q s , the small -x factorization approach is the appropriate 
framework to compute inclusive hadron production in pH collisions, where the gluon saturation 
effects can be systematically included. Therefore, this process pH —> h + X provides a unique 
opportunity to investigate the interplay between two important QCD dynamical effects (small-x 
resummation and the collinear factorization calculation) in high energy hadronic reactions. 

In our calculations, by including the kinematical constraint and the newly added NLO correction 
terms, we are able to improve the results of SOLO and achieve excellent agreement with the 
forward rapidity RHIC data in dAu collisions. Furthermore, we have significantly improved the 
numerical accuracy of the SOLO package, which allows us to compute forward rapidity observables 
at LHC energy with small uncertainties and obtain remarkable agreement with the forward rapidity 






16 


ATLAS data at ^/T/vjv = 5.02 TeV in the relatively low-pj_ region. These results could be additional 
compelling evidence for the observation of the onset of saturation effects at the LHC. 

Our results provide a benchmark framework for the small-a; saturation calculation for high 
energy processes in p A collisions at the next-to-leading order. Recent theoretical developments 
have revolutionized the test of the saturation physics from the qualitative level to the quantitative 
level. With more and more experimental data available from the LHC, we will be able to tell 
whether and when gluon saturation has an effect at extremely small x and large nucleus mass 
numbers. The theoretical advances in computing these processes at the next-to-leading order 
will be crucial to help identify the signature of gluon saturation phenomena. We expect more 
developments along the line discussed in this paper. 

Phenomenologically, the SOLO program has been developed from a single-purpose program, 
for specific formulas under RHIC conditions only, to a more general-purpose program that can 
easily be adapted to different expressions and produce results for both RHIC and the LHC. During 
the development process, we have found an efficient method to perform the numeric computations 
involve. This will become useful and applicable to other interesting processes. 
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Appendix A: The evaluation of some integrals 


First of all, all the Fourier transform formula used in this paper can be derived with the following 
two identities for u± > 0 together with analytical continuation 


roo 2 i+ a rn + —l 

G(a) ee / d( ± ii + “ J o((±»±) = ~ , + 2' 
Jo U I I -n 


roo 

/ dl±l±Jo(l±u±)ln n l± = 

Jo 


d n G(a) 


da" 


(Al) 

(A2) 


a —>-0 


Let us now provide some essential details for deriving L\(k±), which comes from the following 
integral 


a q N r 


IT- 


I _ y±) 


(x_L - b_ l ) 2 


ln Aq(x _l ~ b±) 2 _ (x_L - 6j_) • (y± - b±) ^ kj_\x± - b±\\y± - 6_l| 


= lim 

p—0 27T 


OsNcS± f dV - ik± . r± 


a 


(2?r) 


,N C S ± f d 2 r_ 


ik 


47r 2 


2vr 


' r± S(r±) ( 


(x ± - b±) 2 (y± - b_ l) 2 
S(r ± ) 

k 2 ,r 2 , \ 2 


c o 

i k^\ 2 f°° dl± kJ 

‘v) A 17 ii Jo(,±r - L) 
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where we have used p as an infrared cutoff to compute the above integration. It is important to 
notice that the above result is independent of regularization scheme, since the whole expression 
is finite. Furthermore, by using the following trick (see also Ref. EH ), we can evaluate L\(k±) 
analytically. Let us define 


1 ( 13 ) = / dr_ L r J L +d J 0 (A:_ L r_ L )exp (- 


Q 2 s 


2 1 +/ 3 r[i + f]L[-i-§,-^] 


Q 


2+/3 


(A4) 


where L [—1 — j, x] is defined as the Multivariate Laguerre Polynomial. It is then trivial to find 
that 


L q i(k ± ) = -4 


/"(O) + 2 In —/'(O) + 
co 



(A5) 


which gives L q i(k±) found above. It is useful to note that L[—l,x\ = e x and its first derivative on 
its first argument L^ 1,0 ^[—1, x] = — [7 e + r[0, x] + lnx] e x . 

Now let us further evaluate L\(k±) in the limit of k± -+ 00 . We can rewrite the integral in 
question as follows 



It is straightforward to show that the leading power expansion of the above integral gives 



Appendix B: Some important Fourier transforms 


Although the original derivation mm of the NLO correction for hadron productions in pA 
collisions are completed in the coordinate space, in order to achieve better numerical accuracy, we 
have used Fourier transforms to convert most of the NLO corrections into momentum space in 
Ref. [IS] when the first version of the SOLO package was developed (see (44| for more details of 
the implementation). We only evaluate a couple of NLO terms, for example in coordinate 

space, since the evaluation can be done pretty accurately even in the coordinate space for forward 
rapidity kinematical region at RHIC and the LHC. 

However, when we try to compare to LHC data at middle rapidity, those terms which are left 
in the coordinate space, suddenly give huge numerical uncertainty. The typical integration which 
poses a challenge to numerical integrations is 

/ <bi) 


In the LHC middle rapidity kinematical region, the allowed region of k± is quite large. Note that 
k± = —, which is the transverse momentum of the produced parton, can be much larger than 
p _|_ which is the transverse momentum of the measured hadron. Although the expression is well- 
defined analytically, the above integration oscillates too fast when k± is large, therefore it causes 
a lot of numerical uncertainty. Luckily, we manage to find a way to convert the above integration 
into momentum space which is much more stable. Using the identity 


d 2 Xj 


In 


2 9 1 

x ± p z 


0 -ik ± -x± _ 


l\h - 2 ^ <2 W> 


k L k 
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we find 


d 2 xj 


(27T 


9 5 1 11 (x_l ) In „ 


c o -ifc.-x, = I f 

vr./ Z 2 


^(*i + i±)- JoC-UW'x) 


(B3) 


It is straightforward to test the above identity and find the momentum space expression stable and 
finite. 

In addition, for the same reason, the new L q \ contribution that we found in this paper, as shown 
in Eq. (15), is even more unstable in the coordinate space at the LHC kinematical region. Again, 
we manage to find a nice trick to convert that expression into momentum space. Using dimensional 
regularization in the MS scheme, we can find the following results 
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(B5) 


The derivation of the first expression can be found in Eq. (A4) of Ref. |53j . while the second 
expression can be computed directly. This trick is inspired by the computations of the Sudakov 
factors in the saturation formalism. Subtracting Eq. (B5) from Eq. (B4) and taking e —> 0 gives 
the identity 


, k 2 R 2 \ 2 
m-r) = 8* 


d 2 Z I 


2 r 
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ln 


(27T ) 2 l\ l 2 ± L 


8(k ± - l ± ) - e~ il± ' R± 
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At the end of the day, one can easily find 


d 2 rj 
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ln ■ 
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7T 


d 2 /1 , k '| 


l 2 


In [d(k ± - l±)F{k±) - F(k± + l ± )] . (B7) 


Up to this point, we have transformed all the NLO corrections into momentum space in SOLO 
package, which is relatively more stable numerically at both RHIC and the LHC kinematical region. 
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